/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		subdivide.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-03-22 02:28
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-07-02 12:37
*
*   Description:
*
=============================================================================*/
#ifndef SUB_DIVIDE_H
#define SUB_DIVIDE_H

#include "../../xfem/crack.hpp"
#include "../shapes.hpp"

namespace geoxfem {
namespace algorithm {

template <typename Geometry1, typename Geometry2>
struct subdivide_impl {
  static bool apply(Geometry1& geometry1, Geometry2& geometry2) {
    return SubDivide(geometry2, geometry1);
  }
};
template <typename Geometry1, typename Geometry2>
struct subdivide_without_child_impl {
  static bool apply(Geometry1& geometry1, Geometry2& geometry2) {
    return SubDivide(geometry2, geometry1);
  }
};

template <typename Geometry1, typename Geometry2>
inline bool SubDivide(Geometry1& geometry1, Geometry2& geometry2) {
  return subdivide_impl<Geometry1, Geometry2>::apply(geometry1, geometry2);
}

template <typename DimensionType, typename degree, size_t dim>
struct subdivide_without_child_impl</*{{{*/
                                    geometry::Point<dim>,
                                    geometry::Segment<DimensionType, degree>> {
  using Segment    = geometry::Segment<DimensionType, degree>;
  using SubShape   = geometry::ShapeInterface<DimensionType, degree>;
  using SubSegment = geometry::Segment<DimensionType, degree>;
  using point_type = typename traits::point_type<DimensionType>::type;

  static bool apply(point_type& pt, Segment& seg) {
    /* sub divide the segment */
    point_type p1 = seg.getHeadPoint();
    if(boost::geometry::comparable_distance(p1, pt) > 1e-5) {
      std::vector sub_plist{{p1, pt.deepClone()}};
      SubShape*   sub = new SubSegment{sub_plist};
      seg.addChild(sub);
    }
    point_type p2 = seg.getTailPoint();
    if(boost::geometry::comparable_distance(p2, pt) > 1e-5) {
      std::vector sub_plist{{pt.deepClone(), p2}};
      SubShape*   sub = new SubSegment{sub_plist};
      seg.addChild(sub);
    }
    return true;
  }
}; /*}}}*/

template <typename DimensionType, typename degree1, typename degree2>
struct subdivide_without_child_impl /*{{{*/
  <geometry::Segment<DimensionType, degree1>,
   geometry::PolygonInterface<DimensionType, degree2>> {
  using Segment     = geometry::Segment<DimensionType, degree1>;
  using Polygon     = geometry::PolygonInterface<DimensionType, degree2>;
  using SubShape    = geometry::ShapeInterface<DimensionType, mpl1>;
  using SubSegment  = geometry::Segment<DimensionType, mpl1>;
  using SubTriangle = geometry::Triangle<DimensionType, mpl1>;
  using point_type  = typename traits::point_type<DimensionType>::type;

  static bool apply(Segment& seg, Polygon& poly) {

    point_type seg_p1 = seg.getHeadPoint();
    point_type seg_p2 = seg.getTailPoint();

    std::vector<point_type> intersections;
    boost::geometry::intersection(seg, poly, intersections);

    /* option 0. tip inside the cell */
    if(intersections.size() == 1) {
      point_type* p_tip;
      if(boost::geometry::within(seg_p1, poly))
        p_tip = &seg_p1;
      if(boost::geometry::within(seg_p2, poly))
        p_tip = &seg_p2;

      std::vector<point_type> plist = poly;
      if(plist.back() == plist.front())
        plist.pop_back();
      plist.emplace_back(intersections.front());
      std::sort(plist.begin(), plist.end(),
                algorithm::clock_wise<point_type>(plist));
      plist.emplace_back(plist.front());
      for(auto it1 = plist.begin(); it1 != std::prev(plist.end()); ++it1) {
        auto      it2 = std::next(it1);
        SubShape* sub = new SubTriangle{{*it1, *it2, *p_tip}};
        poly.addChild(sub);
      }
    }
    /* option 1. cut through */
    else {
      GEOXFEM_ASSERT(intersections.size() == 2);
      point_type              ip1   = intersections.front();
      point_type              ip2   = intersections.back();
      std::vector<point_type> plist = poly, sub_plist1, sub_plist2;
      if(plist.back() == plist.front())
        plist.pop_back();
      std::copy_if(plist.begin(), plist.end(),
                   std::inserter(sub_plist1, std::begin(sub_plist1)),
                   [&ip1, &ip2](const point_type& p) {
                     return if_at_left_wing(ip1, ip2, p);
                   });
      std::set_difference(plist.begin(), plist.end(), sub_plist1.begin(),
                          sub_plist1.end(), std::back_inserter(sub_plist2));

      sub_plist1.emplace_back(ip1);
      sub_plist1.emplace_back(ip2);
      point_type pc1 = geometry::getCentroid(sub_plist1);
      std::sort(sub_plist1.begin(), sub_plist1.end(),
                algorithm::clock_wise<point_type>(pc1));
      sub_plist1.emplace_back(sub_plist1.front());
      for(auto it1 = sub_plist1.begin(); it1 != std::prev(sub_plist1.end());
          ++it1) {
        auto      it2 = std::next(it1);
        SubShape* sub = new SubTriangle{{*it1, *it2, pc1}};
        poly.addChild(sub);
      }

      sub_plist2.emplace_back(ip1);
      sub_plist2.emplace_back(ip2);
      point_type pc2 = geometry::getCentroid(sub_plist2);
      std::sort(sub_plist2.begin(), sub_plist2.end(),
                algorithm::clock_wise<point_type>(pc2));
      sub_plist2.emplace_back(sub_plist2.front());
      for(auto it1 = sub_plist2.begin(); it1 != std::prev(sub_plist2.end());
          ++it1) {
        auto      it2 = std::next(it1);
        SubShape* sub = new SubTriangle{{*it1, *it2, pc2}};
        poly.addChild(sub);
      }
    }
    return true;
  }
}; /*}}}*/

template <typename DimensionType, typename degree, size_t dim>
struct subdivide_impl /*{{{*/
  <geometry::Point<dim>, geometry::Segment<DimensionType, degree>> {

  using Segment    = geometry::Segment<DimensionType, degree>;
  using SubShape   = geometry::ShapeInterface<DimensionType, mpl1>;
  using SubSegment = geometry::Segment<DimensionType, mpl1>;
  using point_type = typename traits::point_type<DimensionType>::type;

  static bool apply(point_type& pt, Segment& seg) {

    if(boost::geometry::comparable_distance(pt, seg)
       > boost::geometry::length(seg) * 1e-7)
      return false;

    /* find the sub child containing the point */
    if(seg.hasChild()) {
      auto& children = seg.getChildren();
      auto  it       = std::min_element(
        children.begin(), children.end(),
        [&pt](SubShape* shp1, SubShape* shp2) {
          return boost::geometry::comparable_distance(pt, shp1->toSegment())
                 < boost::geometry::comparable_distance(pt, shp2->toSegment());
        });
      SubShape* child = *it;
      subdivide_without_child_impl<point_type, Segment>::apply(
        pt, child->toSegment());
      seg.deleteChild(child);
    }
    else
      subdivide_without_child_impl<point_type, Segment>::apply(pt, seg);
    return true;
  }
}; /*}}}*/

template <typename DimensionType, typename degree1, typename degree2>
struct subdivide_impl /*{{{*/
  <geometry::Vertex<DimensionType, degree1>,
   geometry::Segment<DimensionType, degree2>> {

  using Vertex     = geometry::Vertex<DimensionType, degree1>;
  using Segment    = geometry::Segment<DimensionType, degree2>;
  using point_type = typename traits::point_type<DimensionType>::type;

  static bool apply(Vertex& vertex, Segment& seg) {
    return SubDivide(vertex.getHeadPoint(), seg);
  }
}; /*}}}*/

template <typename DimensionType, typename degree1, typename degree2>
struct subdivide_impl /*{{{*/
  <geometry::Segment<DimensionType, degree1>,
   geometry::Segment<DimensionType, degree2>> {

  using Segment1   = geometry::Segment<DimensionType, degree1>;
  using Segment2   = geometry::Segment<DimensionType, degree2>;
  using point_type = typename traits::point_type<DimensionType>::type;

  static bool apply(Segment1& seg1, Segment2& seg2) {
    std::vector<point_type> intersection;
    boost::geometry::intersection(seg1, seg2, intersection);
    if(intersection.empty())
      return false;

    GEOXFEM_ASSERT(intersection.size() == 1);

    point_type& ip = intersection.front();
    SubDivide(ip, seg1);
    SubDivide(ip, seg2);

    return true;
  }
}; /*}}}*/

template <typename DimensionType, typename degree>
struct subdivide_impl /*{{{*/
  <geometry::Segment<DimensionType, degree>,
   geometry::PolygonInterface<DimensionType, degree>> {

  using Shape      = geometry::ShapeInterface<DimensionType, degree>;
  using SubShape   = geometry::ShapeInterface<DimensionType, mpl1>;
  using Segment    = geometry::Segment<DimensionType, degree>;
  using Polygon    = geometry::PolygonInterface<DimensionType, degree>;
  using point_type = typename traits::point_type<DimensionType>::type;

  static bool apply(Segment& seg, Polygon& poly) {
    bool if_divided = boost::geometry::covered_by(seg, poly)
                      || boost::geometry::crosses(seg, poly);
    if(!if_divided)
      return false;

    /* find the sub child cutted by the segment */
    if(poly.hasChild()) {
      auto& children = poly.getChildren();
      auto  it       = std::find_if(
        children.begin(), children.end(), [&seg](const SubShape& shp) {
          return boost::geometry::covered_by(seg, shp.toPolygon())
                 || boost::geometry::crosses(seg, shp.toPolygon());
        });
      SubShape* child = *it;
      subdivide_without_child_impl<Segment, SubShape>::apply(seg, child);
      poly.deleteChild(child);
    }
    else
      subdivide_without_child_impl<Segment, SubShape>::apply(seg, poly);

    return true;
  }
}; /*}}}*/

template <typename DimensionType, typename degree>
struct subdivide_impl</*{{{*/
                      geometry::Segment<DimensionType, degree>,
                      geometry::Mesh<DimensionType, degree>> {
  using Mesh    = geometry::Mesh<DimensionType, degree>;
  using Shape   = geometry::ShapeInterface<DimensionType, degree>;
  using Segment = geometry::Segment<DimensionType, degree>;

  static bool apply(Segment& crack_seg, Mesh& mesh) {

    std::set<Shape*> cell_list;

    for(auto& seg_in_mesh : shapesIn(mesh, 1))
      if(SubDivide(crack_seg, seg_in_mesh->toSegment())) {
        seg_in_mesh->setXfem(true);
        for(auto cell : neighborsIn(*seg_in_mesh, 2)) {
          cell->setXfem(true);
          cell_list.emplace(cell);
          for(auto vertex : neighborsIn(*cell, 0))
            vertex->setXfem(true);
        }
      }
    model::buildGaussPoints(crack_seg);

    for(auto& cell : cell_list) {
      bool if_divided = SubDivide(crack_seg, cell->toPolygon());
      if(if_divided) {
        model::buildGaussPoints(*cell);
        cell->updateGaussXfem(crack_seg);
      }
    }
    return true;
  }
}; /*}}}*/

template <typename DimensionType, typename degree>
struct subdivide_impl</*{{{*/
                      xfem::Crack<DimensionType, degree>,
                      geometry::Mesh<DimensionType, degree>> {
  using Crack   = xfem::Crack<DimensionType, degree>;
  using Mesh    = geometry::Mesh<DimensionType, degree>;
  using Segment = geometry::Segment<DimensionType, degree>;

  static bool apply(Crack& crack, Mesh& mesh) {

    for(auto& crack_seg : shapesIn(crack, 1))
      SubDivide(crack_seg->toSegment(), mesh);

    return true;
  }
}; /*}}}*/

}  // namespace algorithm
}  // namespace geoxfem

#endif /* SUB_DIVIDE */
